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ABSTRACT 

Directly imaging extrasolar terrestrial planets necessarily means contending 
with the astrophysical noise of exozodiacal dust and the resonant structures cre- 
ated by these planets in exozodiacal clouds. Using a custom tailored hybrid 
symplectic integrator we have constructed 120 models of resonant structures cre- 
ated by exo-Earths and super-Earths on circular orbits interacting with colli- 
sionlcss steady-state dust clouds around a Sun-like star. Our models include 
enough particles to overcome the limitations of previous simulations that were 
often dominated by a handful of long-lived particles, allowing us to quantitatively 
study the contrast of the resulting ring structures. We found that in the case 
of a planet on a circular orbit, for a given star and dust source distribution, the 
morphology and contrast of the resonant structures depend on only two param- 
eters: planet mass and v /o^//3, where a p is the planet’s semi-major axis and (3 
is the ratio of radiation pressure force to gravitational force on a grain. We con- 
structed multiple-grain-size models of 25,000 particles each and showed that in a 
collisionless cloud, a Dohnanyi crushing law yields a resonant ring whose optical 
depth is dominated by the largest grains in the distribution, not the smallest. 

We used these models to estimate the mass of the lowest-mass planet that can 
be detected through observations of a resonant ring for a variety of assumptions 
about the dust cloud and the planet’s orbit. Our simulations suggest that planets 
with mass as small as a few times Mar’s mass may produce detectable signatures 
in debris disks at a p > 10 AU. 
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1. Introduction 


A number of proposed experiments like the Terrestrial Planet Finder (TPF) aim to 
directly image the scattered and emitted light from extrasolar planets (Lawson & Traub 
2006). These experiments will also excel at detecting exozodiacal dust, circumstellar dust 
analogous to zodiacal dust in our solar system (e.g. Agol 2007; Beckwith 2007). Zodiacal 
dust in the solar system consists of ~ 1 — 100 /irn dust grains released through asteroidal 
collisions and the outgassing of comets (e.g. Schramm et al. 1989). This dust forms the 
zodiacal cloud, extending from the solar corona (e.g. Mann et al. 2000) to beyond Jupiter 
(e.g. Kruger et al. 1999). 

Our zodiacal cloud exhibits several structures interpreted as dynamical signatures of 
planets (Dermott et al. 1985, 1994; Reach et al. 1995). Several dusty disks around nearby 
main-sequence stars show similar structures (e.g. Greaves et al. 1998; Wilner et al. 2002; 
Kalas et al. 2005). This trend suggests that exozodiacal clouds may be full of rings, clumps 
and other asymmetries caused by planets and other phenomena. 

This situation raises some important questions. Will the structures in exozodiacal clouds 
be harmful astrophysical noise for direct imaging of extrasolar planets (Beichman 1996; 
Beichman et al. 1999)? Or can the dynamical signatures of planets in these clouds help us 
find otherwise undetectable planets (e.g. Kuchner & Holman 2003)? 

Several studies have examined the geometry of resonant signatures of planets in debris 
disks (e.g. Kuchner & Holman 2003; Reche et al. 2008). However, most simulations cannot 
quantitatively study the contrast in these structures: how bright they are relative to the 
background cloud. We need to model the contrast of the structures in exozodiacal clouds to 
understand their roles as astrophysical noise and as signposts of hidden planets. However, 
accurately simulating the contrast of these structures demands computational resources that 
have only recently become available (e.g. Dcller & Maddison 2005). 

In this paper we examine the contrast of resonant structures induced by planets in 
steady-state exozodiacal clouds and the detectability of these structures via direct imaging. 
We simulate high-fidelity images of collisionless exozodiacal clouds containing a terrestrial- 
mass planet — an exo-Earth or super-Earth. By using roughly an order of magnitude more 
particles than most previous simulations, we overcome the Poisson noise associated with 
constructing histograms of the column density and populating the external mean motion 
resonances (MMRs) of planets. We use our simulations to estimate the minimum planet 
mass that can be indirectly detected via observations of these structures as a function of 
the planet semi-major axis and dominant grain size under the assumption of circular planet 
orbits. Our models apply to exozodiacal clouds less than a few hundred times the optical 
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depth of the solar zodiacal cloud, clouds for which the collision time is shorter than the 
Poynting-Robertson (PR) time for typical grains. 

Section 2 of this paper describes our numerical techniques. We present a synthetic 
catalog of resonant debris disk structures in Section 3. We describe our multiple-particle- 
size cloud models and discuss their detectability in Section 4. In Section 5, we discuss the 
limitations of our simulations; we summarize our conclusions in Section 6. 


2. Numerical Method 

Dust grains in the inner solar system are primarily released from parent bodies via colli- 
sions or outgassing. Radiation pressure ejects the smallest particles from the solar system in 
a dynamical time while the larger particles slowly spiral inward due to PR drag (Robertson 
1937; Burns et al. 1979). During their spiral toward the Sun, particles may become tem- 
porarily trapped in the MMRs of planets, extending their lifetimes by a factor of a few to 
ten (Jackson & Zook 1989). This trapping locally enhances the particle density, creating 
structures within the zodiacal cloud, which have been described as circumsolar rings, bands, 
and clumps (e.g. Kelsall et al. 1998). 

To model these types of structures in exozodiacal clouds we numerically integrated the 
equation of motion of dust particles. The equation of motion for a perfectly absorbing 
particle orbiting a star of mass m* is given to first order in v/c by Robertson (1937): 


d 2 r 
dt 2 


(l + sw)/R?m*,. A , , 

2~(1 - P)r 5 “ [rr + v , 




C 




( 1 ) 


where r and v are the heliocentric position and velocity of the particle and sw is the ratio 
of solar wind drag to PR drag. We assume a value for sw of 0.35 (Gustafson 1994). For 
perfectly absorbing spherical particles in the vicinity of the Sun, f3 ~ 0.57 /ps, where p is the 
mass density of the particle in g cm -3 and s is the radius in //m. 


2.1. A Customized Hybrid Symplectic Integrator 

We implemented a customized hybrid symplectic integrator to perform our numerical 
integrations. Chambers (1999), hereby referred to as C99, introduced hybrid symplectic 
integration as a method for dealing with close encounters in an efficient n-body code. Sym- 
plectic integrators rely on splitting the Hamiltonian into two easily integrable portions — a 
dominant term, Hp, and a smaller perturbative term, Hp. However, in the n-body problem, 
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Hp may exceed Hp, during close encounters. Hybrid symplectic integrators overcome this 
problem by effectively switching from a symplectic integrator to an alternate integrator (e.g. 
Bulirsch-Stoer). 

The hybrid method reduces the perturbative term of the Hamiltonian, iip, by a factor 
K(rij), where r l} is the distance between the two bodies in question, to ensure that the 
perturbative term remains relatively small. The integrator includes the remaining portion of 
the perturbative term, . Hp tJ [\ — K (ry,-)], in the dominant term which is then integrated 
using a method of choice. The “changeover function,” K(r t j), is a smooth function that 
varies from 0 for nj < r crit to unity for r tJ > r crit . 

Using a hybrid integrator requires choosing a changeover function and a value for r cr it. 
We use the same changeover function as C99. We assign a different value of r crit ti to each 
body, calculated as the larger of 3i?n,i and Try, where i? H ,i and ry are the Hill radius and 
velocity of the i th body, respectively, and r is the time step of the integrator. We then 
calculate the critical distance for a pair of bodies as r cr i t ,p- = r cr it ; * + r cr i t j. 

Our integrator also incorporates the effects of radiation pressure, PR drag, and solar 
wind drag. We implement radiation pressure as a correction to the effective stellar mass (cf. 
Eq. 1) and treat the drag effects as an additional term in Hp, in much the same way as 
Moro-Martin & Malhotra (2002), hereby referred to as MMM02. We also use democratic 
heliocentric (DH) coordinates, composed of the barycentric momenta and heliocentric posi- 
tions, because of their relative ease of implementation. This choice introduces an additional 
perturbative term to the Hamiltonian due to the motion of the star with respect to the 
barycenter (Duncan et al. 1998). 


2.2. Comparison of Integrator with Previous Results 

We checked our integrator using a variety of standard tests. We checked the energy and 
Jacobi constant conservation with the drag terms turned off and examined the evolution of 
dust particles’ orbital elements under our implementation of drag effects. We also compared 
our hybrid integrator to a Bulirsch-Stoer integrator by examining the path of an individual 
test particle during a close encounter and by examining the statistics of a cloud of particles 
in a collisionless disk containing a planet. 

We tested energy conservation in our integration code by integrating the orbits of the 
four outer planets and the Sun for 3 x 10 5 years using a time step of 0.15 years. The energy 
error was bounded with a mean value of AE/E « 3 x 10~ 9 . 
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Duncan et al. (1998), which we will refer to as DLL98, tested the relative conservation 
of energy in their symplectic integrator as a function of planet perihelion distance. We 
replicated their tests using our code. Figure 1 shows the relative energy error in an integration 
of the orbit of Jupiter for 3 x 10 5 years using a time step of 0.15 years. We initially placed 
Jupiter at aphelion. Figure 1 also shows the results of integrating the orbits of Jupiter 
and Saturn under the same conditions. With the DH method, the perturbative solar term 
increases as the perihelion distance of the planet decreases, causing the fractional energy 
error to increase similarly. The fractional energy errors shown in Figure 1 agree with those 
obtained by DLL98. 

We checked the conservation of the Jacobi constant by integrating particles in the Sun- 
Neptune system. We found results consistent with those of MMM02. Particles that did not 
undergo close encounters conserved the Jacobi constant at the level of rs_/ 10“ 8 to 10“ 7 . 

To test our implementation of PR drag, we replicated a test performed by MMM02. 
We integrated the orbit of a particle with f3 = 0.2 and sw = 0.35 in the presence of the 
Sun. Figure 2 shows the semi-major axis and eccentricity as functions of time. These results 
match the results of MMM02 and agree with the analytic solution (Wyatt & Whipple 1950). 

We tested the performance of our hybrid scheme by integrating the orbit of comet 
P/Oterma in a close encounter with Jupiter, which has been done previously by Michel & 
Valsecchi (1996) and C99. The initial conditions for both bodies can be found in Table 3 of 
Michel & Valsecchi (1996). Figure 3 shows the path of comet P/Oterma for several values 
of integration time step r as seen in the frame co-rotating with Jupiter, which is located 
at the origin. These results are similar to those obtained by C99. Our code shows a minor 
improvement over the other codes, most noticeable in the r = 100 days case, that is likely 
only clue to differences in the calculation of the changeover distance. C99 explicitly sets 
v crit = 3 J?h f° r this test; we used our prescription for r crit as described in Section 2.1. 


2.3. Test Simulations of a Steady-State Exozodiacal Cloud 

We directly compared simulations of resonant structures made with our hybrid integra- 
tor to simulations made with a Bulirsch-Stoer integrator. During the integrations we recorded 
the coordinates of each particle in a 2-D histogram at regular intervals. This histogram mod- 
els the surface density distribution in a steady-state cloud. Since we only modeled planets on 
circular orbits, we simply recorded the coordinates in the frame co-rotating with the planet. 
This technique has been widely used by dust cloud modelers (Dermott et al. 1994; Liou & 
Zook 1999; Moro-Martm & Malhotra 2002; Wilner et al. 2002; Dellcr & Maddison 2005). 
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Figure 4 shows two histograms, one for each integrator, for simulations of 1,000 particles 
each in the presence of the Sun/Earth system. We used a histogram bin width of 0.0175 AU. 
For these simulations we chose (3 = 0.02 and initially released the particles with semi-major 
axis, G^ust, distributed uniformly between 3 and 5 AU, eccentricity, e, uniformly distributed 
from 0.0 to 0.1, and inclination, i, uniformly distributed between 0° and 6°. We used a 
symplectic time step of 0.02 years and recorded the particle locations every 250 years. 

Except for a small number of pixels, the middle panel of the figure (simulation using the 
hybrid integrator) looks qualitatively very similar to the left-most panel (simulation using 
the Bulirsch-Stoer integrator). The right panel of Figure 4 shows the difference of these two 
images divided by the ^Jn Poisson noise expected for each pixel where n is the number of 
particles in the pixel. This figure demonstrates that the differences between the two models 
are nearly consistent with the Poisson noise of the histograms. The two integrators resulted 
in histograms with minor structural differences, but the hybrid symplectic integrator runs a 
few times faster. 

Besides pixcl-to-pixel Poisson noise in the histogram, this method is also sensitive to 
noise in the population of MMRs. MMM02 showed that the population of the dominant 
resonances varied by a factor of ~ 3 among sets of 100-particle simulations of Kuiper Belt 
dust interacting with Neptune. This noise probably causes the differences between the two 
simulations shown in Figure 4 beyond those attributable to pixcl-to-pixel Poisson noise. 
Although simulations of 100 particles may acquaint us with the generic geometry of debris 
disk structures, we cannot use them to predict ring contrasts; to model the contrast in a 
resonant clond feature we must include enough particles to accurately populate the MMRs. 

We solved this problem by using more particles. We used the 420-processor Thunderhead 
cluster at NASA Goddard Space Flight Center to perform simulations of 5,000 particles each. 
Figure 5 shows the population of MMRs for three independent 5,000-particle simulations of 
the Sun and four outer planets using the same initial conditions as MMM02. Simulating 
5,000 particles reduced the difference between MMR populations for the three simulations 
to less than 7% for the dominant 2:1 and 3:2 MMRs, allowing us to synthesize high-fidelity 
images and quantitatively study the resonant ring structures. 


3. Simulations & Results 

3.1. Cataloging Debris Disk Structure 

To explore the range of different types of structures formed by terrestrial- mass planets, 
we performed 120 simulations of dust interacting with single planets on circular orbits. The 



simulations used 5,000 particles each and covered six values of planet mass, M p (0.1, 0.25, 
0.5, 1.0, 2.0, and 5.0 M e ), four values of planet semi-major axis, a v (1, 3, 6, and 10 AU), and 
five values of ft (0.0023, 0.0073, 0.023, 0.073, and 0.23) corresponding to spherical silicate 
particles ranging in radius from ~ 1 — 120/jrn. We released the particles on orbits with 
semi-major axes uniformly distributed between 3.5 and 4.5 times the semi-major axis of the 
planet’s orbit — well outside of the strongest MMRs. We used initial eccentricities uniformly 
distributed between 0 and 0.2, initial inclinations uniformly distributed between 0 and 20°, 
and the longitude of the ascending node, fl, and the argument of pericenter, c o, uniformly 
distributed between 0 and 2ti. We considered planet semi-major axes of 1 to 10 AU because 
typical designs for TPF can detect an exozodiacal cloud with 10 times the optical depth of 
the solar zodiacal cloud over roughly that range of circumstellar radii (Levine et al. 2006). 

We chose these initial conditions to model only dynamically-cold dust, i.e. ed us t 0.2 
and idust < 20°, since this component of a dust cloud is the dominant contributor to resonant 
ring structure. We neglect dynamically hot dust with the idea that it can always be added 
in later as a smooth background (Moran et al. 2004). The asteroid belt probably produces 
much of the solar system’s dynamically cold dust, while comets are thought to contribute a 
more dynamically hot cloud component (e.g. Liou et al. 1995; Ipatov et al. 2007). We treat 
only steady-state dust clouds, assuming dust is continually replenished, and ignore transient 
collisional events. 

Figure 6 shows some examples of the histograms from our simulations, which reveal a 
wide range of trapping behavior. Some histograms show no azimuthal or radial structure, 
while others show high contrast rings. All of the patterns are Type I structures as identified 
by Kuchner & Holman (2003). 

Several general trends emerged. The ring contrast increased with increasing planet 
mass. Reducing ft also enhanced trapping, as did increasing the planet’s semi-major axis. 

These last two trends can be explained by comparing the libration time of a given MMR 
to the PR time. The PR time scales as a| ust //3, while the libration time for a given resonance 
scales as a^ t , where adust is the semi-major axis of a dust grain’s orbit. The ratio of these 
quantities yields ^ adust /ft, a parameter that measures the degree to which resonant trapping 
is adiabatic (e.g. Henrard 1982); the trapping becomes more adiabatic and more efficient at 
greater distances from the star, and for larger particles. We discuss this phenomenon further 
in Section 3.3 below. 

In addition to following these trends, all of the simulated ring structures, like those 
shown in Figure 6, share some salient features: 

1. For cases in which even a modest amount of trapping occurs (azimuthally averaged 



contrasts ^ 1.3 : 1), the ring structures exhibit a sharp inner edge at ~ 0.83a p . This 
feature probably appears because the eccentricities of particles trapped in exterior 
MMRs are typically pumped up to a limiting value before a close encounter with the 
planet ejects them from resonance. For a particular MMR, all particles, regardless of 
/3, tend to approach the same limiting eccentricity and accordingly, a similar percen- 
ter distance (Beauge & Ferraz-Mello 1994). The limiting eccentricities are such that 
the limiting pericenter distances are nearly equal for the dominant resonances (e.g. 
the 2:1 and 3:2 resonances have limiting pericenter distances of 0.823a p and 0.827a p , 
respectively), creating the ring structure’s sharp inner edge. 

2. A gap in the ring structure, a local minimum in the surface density, appears around 
the planet. If we define gap width as the FWHM of the minimum in the azimuthal 
surface density profile at r = a p , we find that the gap width is linearly proportional 
to the contrast of the ring, as shown in the left panel of Figure 7. A linear fit to the 
data shown in this figure gives w gap ~ 10° x Caa,ie for Caa.ie > 1-6, where w gap is the 
gap width in degrees and Caa,ie is the azimuthally averaged inner-edge contrast (see 
Section 3.2). 


3. 


The rings show a leading-trailing asymmetry. The trailing side of the ring structure 
is noticeably denser than the leading side, and the structure is rotationally shifted in 
the prograde direction causing the trailing side to be closer to the planet (Dermott et 
al. 1994). To examine the leading-trailing asymmetry caused by a prograde shift of 
the ring structure, we measured the azimuthal offset of the center of the gap described 
above from the planet. The right panel of Figure 7 shows these measured prograde 
shifts of our simulations. Kuchner & Holman (2003) showed for a particular first order 
exterior MMR, 
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where (j ) r ing is the measured prograde shift of the ring structure. While the relationship 
in Equation 2 holds for a single MMR, it does not strictly apply to a given ring 
structure which consists of several well-populated MMRs. The relative populations of 
these MMRs are also functions of M p , a p , and /?. However, this situation seems to 
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preserve a power-law relationship between sin0 ring and (3 (1 — (5) / (M PV /oj/), as shown 
in Equation 3. 

4. The radial width of the ring increases with the contrast of the ring, ranging from a few 
percent of a p to ~ 1.6a p in the highest contrast case. As the trapping probabilities of 
all the MMRs increase, MMRs farther from the planet’s orbit become populated. For 
this reason, the outer-edge of the ring structure differs significantly among simulations. 
The outer-edge can be quite blurry or very well-defined, making the radial width of a 
ring structure difficult to quantify. 

Our catalog of debris disk structures induced by terrestrial-mass planets is publicly avail- 
able online at http://asd.gsfc.nasa.gov/Christopher.Stark/catalog.php. This online catalog 
also contains images synthesized from the density distributions in scattered light and 10 /jrn 
thermal emission assuming blackbody grains. Future studies of resonant ring structures with 
TPF or other experiments can use our catalog to interpret dust cloud patterns in terms of 
planet and dust parameters, assuming the observed image is dominated by a single grain 
size. We envision the following process, inspired by recent papers on disks observed with the 
Hubble Space Telescope (e.g. Clampin et al. 2003; Kalas et al. 2005): 

1. Deproject the image to remove inclination effects. 

2. Remove any smooth backgrounds by a power law fit. 

3. Estimate the dominant grain size in the resonant ring using infrared photometry or 
other methods. 

4. Compare the image of the disk to the online catalog to constrain the planet’s mass and 
location. 


3.2. Ring Contrast 

We considered three different metrics for describing the ring contrast in our simulations: 

Cuax'- The surface density of the ring at its densest point divided by the surface density 
of the background cloud 

CwA.Max- The maximum value of the azimuthally averaged surface density divided by 
the surface density of the background cloud 

Caa,ie: The azimuthally averaged surface density at the inner edge of the ring divided 
by the surface density of the background cloud 



10 


We calculated the above contrast metrics for all 120 simulations. We measured the 
surface density of the background cloud at a circumstellar distance r « 0.8a p . The surface 
density of the background cloud was nearly constant inside and outside of the ring, but did 
exhibit a small local minimum near r cs 0.8a p in a few cases. 

To calculate Cmsx, we must search for the densest pixel, which introduces a bias toward 
pixels that exhibit an extreme amount of Poisson noise. To reduce this noise we averaged the 
surface density over nine pixels centered on the densest point. Using C*AA,Max or Caa,ie, on 
the other hand, automatically averages over the effects of Poisson noise in onr simulations. 

Figure 8 shows two examples of how the contrast, Uaa,ie, depends on planet mass and f3. 
Both plots show a similar behavior with three distinct regions: a no-trapping regime (contrast 
~1), a transitional regime, and a saturation regime (maximum contrast). The saturation 
regime is of particular significance. Our results suggest that within the range of parameters 
investigated, for a given value of (3, all contrasts converge to the same value for large planet 
masses independent of planet semi-major axis, i.e. the contrast becomes “saturated” and 
increasing the planet’s semi-major axis has little effect on the contrast. The right panel in 
Figure 8 illustrates this behavior; all four contrast curves, each of which corresponds to a 
different planet semi-major axis, approach the same value of ~ 7 near M p = 5 M©. Similarly, 
for a given planet mass, contrasts converge to the same value for small /3 independent of a p , 
as shown in the left panel in Figure 8. The morphology of the structure can vary, but the 
contrast of the ring structure is roughly constant in these saturation regimes. 


3.3. Adiabaticity 

As we mentioned above, dividing the PR time by the libration time of a given MMR 
yields a parameter, that indicates the degree to which the resonant trapping is 

adiabatic. We plot the contrast in our simulations as a function of this parameter in Figure 
9. This figure demonstrates that for M p < 5 M© and for a given distribution of parent body 
orbital elements, the ring contrast is a function of only two parameters: planet mass and 

\f®p/ P- 

The morphology of the resonant rings is also, to good approximation, a function of only 
the planet mass and f3. The models shown in the two right panels of Figure 6 illustrate 
this phenomenon; for both models « 137 AU 1//2 . These two models have the same 

morphology to a level consistent with pixcl-to-pixel Poisson noise. Note that for small (3 in 
Equation 3, the prograde shift is approximately a function of (y^//?) for a given planet 


mass. 
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For large values of (3 and M p , the morphology and contrast of the ring structures are not 
simple functions of y / a^/[3. Simulations with large values of f3, but equal values of y/o^//3 
(e.g. a/1 AU/0.073 and a/10 AU/0.23) show morphological differences, including differences 
in prograde shift. Our simulations with M p = 5 M® also show contrast differences among 
rings with equal values of v /cqj//3. 

Wyatt (2003) investigated resonant trapping in MMRs for a system of planetesimals 
exterior to an outward migrating planet on a circular orbit. Wyatt (2003) plotted the 
trapping probability for a single MMR in his model as a function of migration rate and planet 
mass and found it could be well approximated by a function of the form P = [l + (a p /pi) P2 ] _1 , 
where a p is the migration rate and the parameters p± and P2 are power laws in planet mass. 
Our trapping scenario assumes dust migrating inward toward the planet, but the concept is 
similar. Since contrast is closely related to trapping probability, we decided to fit the data 
shown in Figure 9 with a function of the form 


C 


l+Pi 





( 4 ) 


inspired by Wyatt (2003). Each of the three parameters, p t , is a power law in planet mass 
of the form pi = pi^Mp 1 ' 2 . We fit all 120 contrast measurements with this six-parameter 
function for each of the three contrast metrics. The best fits, two of which are shown in 
Figure 9, are: 

C'aa.ie: Pi ~ 4.38 (M p /M®) 0 ' 19 , p 2 « 207 (M p /M e )” L17 , p 3 « 2.05 (M p /M®) an 

C'AA.Max: Pi « 4.54 (M p /M®) 0 ' 17 , p 2 « 205 (Mp/M®)" 1 ’ 17 , p 3 « 1.63 (M p /M®) a19 

CW Pi ~ 6.23 (M p /M e ) a27 , p 2 « 164 (M p /M e )" L09 , p 3 « 1.72 (M p /M e ) a05 

Equation 4, combined with the above values, summarizes our results for all combinations 
of planet mass, planet semi-major axis and (3 we simulated. Figure 9 shows the inner- 
edge contrast, Caa.ie, deviates significantly from the fits for large M p and large y /a^//3. 
The increased trapping efficiency for MMRs with these massive planets likely enhances the 
population of MMRs farther from the planet’s orbit and depletes the inner MMRs that cause 
the sharp inner edge. 
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Multi-Particle-Size Models 


4.1. Composite Simulations 

We used our 120 simulations to produce 20 multiple-particle-size dust cloud models 
by forming weighted sums of the histograms assuming a Dohnanyi distribution of particle 
sizes (Dohnanyi 1969). Each of these composite models effectively utilizes 25,000 particles. 
Exactly how we apply the ideas in Dohnanyi (1969) has profound effects on our composite 
models, so we present here two different kinds of models. 

First, we assembled a composite model in which the particles arc initially released from 
their parent bodies according to a crushing law, and do not undergo any further collisional 
processing as they spiral inward. This scenario models a sparse disk with a belt of dust- 
producing material, like our own zodiacal cloud. The crushing law for asteroid material at 
micron sizes is unknown, so we choose the crushing law used by Dohnanyi (1969): 


dN 

ds 


oc s 


( 5 ) 


where dN is the number of particles with radius s in a bin of width ds, and a = 3.4. 


We calculate the optical depth, r, for our composite models from r = JT WjAjCq, where 
Wi, Ai, and cq are the weighting factor, particle cross-section, and surface number density 
of the i th single-particle simulation, respectively. We assume that the cross-section of each 
particle is Ai oc /3 -2 . The crushing law in Equation 5 implies a weighting factor for the i th 
histogram of wy = /?" _2 A/lj, where A/3, is the width of the /3, bin. For a constant logarithmic 
spacing in (3, like the spacing we used in our simulations, and the Dohnanyi crushing law, 
Wi = A 2 - 4 . 

Larger particles have longer PR times, so in the absence of collisional processing, their 
density is enhanced by a factor of /3 _1 under our assumption of a steady-state cloud model. 
One might expect that this effect must be included in the weighting factor. However, our 
simulations include this effect automatically as long as we keep the frequency with which 
particle locations are recorded constant among all of our simulations. We did, in fact, vary 
the recording frequency with the PR time, but we corrected for the differences in recording 
frequency before summing the histograms. 

Figure 10 shows the optical depth of one of our 20 composite models (M p = 2.0 M®, 
a p = 6.0 AU), together with the optical depths of single-particle-size models using only the 
smallest and largest particle sizes included in the composite model. Although the crushing 
law used by Dohnanyi (1969) favors smaller particles by number, even more than some 
empirical crushing laws (Durda et al. 2007), the optical depth in the composite models is 
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dominated by the largest particles. This situation occurs because the larger particles are 
both longer lived (tpR oc /? -1 ) and more likely to be trapped in MMRs. Hence, the upper 
left panel of the figure closely resembles the lower right panel. 


Next, for the purpose of illustration, we ignored the initial size distribution of dust 
particles and forced the disk to obey a size distribution of 



at a radius of ~ 3a p from the star. This scenario probably doesn’t have a physical inter- 
pretation, but it illustrates an interesting phenomenon: how resonant trapping tends to sort 
particles by size. We enforce the size distribution at one location within the disk, but the 
size distribution will not follow a Dohnanyi distribution elsewhere in the disk. 


The top right panel in Figure 10 shows the optical depth of an example of this kind of 
composite cloud, normalized to a Dohnanyi distribution at ~ 18 AU. Models constructed in 
this fashion are more greatly affected by the smallest grains. Hence, the top right panel does 
not greatly resemble the lower right panel in Figure 10. 


4.2. Semi- Analytic Treatment 

We can further develop these ideas with a simple semi-analytic treatment. For a given 
planet mass and semi-major axis, the contrast function (see Equation 4) becomes C'(s), 
where s is the particle size. We approximate the contrast function in Equation 4 with the 
piecewise function 

' 

1 for s < Si 

C « = | (yf for Si < s < S 2 < 7 > 

k d’large for S -S' 2, 

where Ci arge = 1 + pi is the contrast for the largest particles, m is the logarithmic slope 
of the contrast in the transition regime, and .sq and S 2 are the particle sizes that mark the 
beginning and end of the transition regime, respectively. We fit our contrast data with this 
piecewise function and obtained the following power law estimates assuming silicate grains 
(p ~ 2 g cm -3 ): 

/ M \ 0 ' 19 

Cjarge;AA,IE ~ 1 + 4.38 ( J (8) 
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In the same manner as Section 3.2, we defined the contrast of any ring structure as the 
surface density within the ring, <7 r i ng , divided by the background surface density, ctbg- The 
contrast in optical depth of a cloud containing several components of various sized particles, 
labeled with the index i, is 


(Cr) = 


yy &bg, iAi 


(12) 


For a collisionless cloud with a continuous distribution of grain sizes, the contrast in optical 
depth of the composite cloud is given by 


s s - a C(s ) ds 


(Cr) = 


s 3 ~ a ds 


(13) 


where we have explicitly included the particle cross section (A(s) oc s 2 ) and background 
surface density (obg(s) oc s 1_q ). For a collisionless cloud, the background surface density is 
enhanced by a factor of s due to the PR time scaling as s (see Section 4.1), and a factor of 
s~ a , which describes the assumed crushing law. 

Using Equations 7, we can now integrate Equation 13 directly. Assuming s min < si < 
s 2 < W and (s min /s max ) |4_a| <C 1 when 4, we find 
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For cases in which the maximum particle size in a disk is less than s 2 (see Equation 11), 
simply replace all instances of s 2 in Equations 14 with s max . Similarly, for cases in which the 
minimum particle size in a disk is greater than Si (see Equation 10), replace all instances of 
Si with s min . 


Equations 14, together with Equations 8-11, give analytic expressions for optical depth 
contrast in terms of M p , a p , .s min , and s max . Although Equations 14 address all possible 
scenarios, the most plausible scenarios have crushing laws with a < 4 (e.g. Durda et al. 
2007). With this assumption, Equations 14 combined with Equations 8-11 gives 
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If, as in our first composite model in Section 4.1, we assume that the particles are 
released from their parent bodies in accordance with the Dohnanyi (1969) crushing law and 
then spiral inward without colliding, a = 3.4. In this case, Equations 14 give a contrast in 
optical depth of (C T ) ~ Ci ai .g e in the limit s max 3> sy • This result confirms our numerical 
results for our first composite model, shown in the upper left panel of Figure 10; the contrast 
in optical depth is dominated by the large particles. 


For our second composite model, we forced the background density to obey a Dohnanyi 
distribution at ~ 3a p , i.e. ctbg(s) oc s -35 , so that a = 4.5. This technique essentially 
removes the factor of s in the background surface density that results from the PR time 
scaling as s. For the composite cloud shown in Figure 10, M p = 2.0 M e , and a p = 6.0 
AU, for which m ~ 0.68, sy « 1.9 /irn, s 2 ~ 24 /mi, and Ci argei AA,iE ~ 6. We let each 
simulated particle size represent a range of particle sizes using the midpoint method, which 
gives s m i n ~ 0.7 /mi. With these values, Equations 14 give a contrast in optical depth of 
(CV,aa,ie) ~ 2.4, in agreement with the measured contrast in the top right panel of Figure 
10 . 


Figure 11 illustrates in general how a distribution of particle sizes affects the contrast of 
a ring structure. This figure compares the contrast of a collisionless multi-particle-size cloud 
(Equations 14) to that of a single-particle-size cloud as a function of / [4 m - m assuming a 



16 


Dohnanyi (1969) crushing law. Both kinds of clonds have the same contrast in the adiabatic 
limit (large but the contribution of the smaller grains reduces the contrast else- 

where, effectively broadening the transition between the no-trapping regime and saturation 
regime. Crushing laws with a < 3.4 result in contrast curves that more closely resemble the 
single-particle-size contrast curves shown in Figure 11. 

In a real zodiacal cloud, collisions affect the distribution of grains, even far from the 
source of the grains. Our composite dust cloud models do not include collisions and become 
unreliable for particles with collisional times less than their PR times. Our composite models 
also lack the structural results of collisional effects, such as the loss of particles as a function 
of circumstellar distance (Wyatt 2005) and any potential morphological effects in the ring 
structure. 

More sophisticated models may be required to investigate these phenomena. However, 
since dust produced according to a Dohnanyi (1969) crushing law or a Durda et al. (2007) 
crushing law yields a cloud dominated by the largest grains, as we showed above, we hypoth- 
esize that resonant rings in exozodiacal clouds may often be dominated by a single particle 
size whose PR time is roughly equal to the its collisional time. 


4.3. Ring Detectability 

The detectability of a resonant ring structure depends on many factors specific to the 
telescope being used and the observing conditions. We address this complicated issue by 
imposing one simplifying assumption: a minimum detectable optical depth ring contrast of 
1.5. This assumption likely underestimates the sensitivity of a TPF-likc mission to rings in 
exozodiacal clouds analogous to the solar zodiacal cloud. In such a cloud, a ring 0.4 AU 
wide located at 1 AU from the star has ~ 15 times the total flux of an Earth-like planet at 1 
AU, even for a contrast of unity. Our assumption, conservative on the basis of photon noise 
alone, allows for the possibility of unknown systematic noise that could hinder the detection 
of extended structures. 

Figure 12 shows the minimum detectable planet mass as a function of semi-major axis 
and maximum dust particle size based on Equations 15 and a Dohnanyi (1969) crushing law. 
The masses and semi-major axes of Earth, Mars, and the planet OGLE-2005-BLG-390Lb, 
detected by the microlensing technique (Beaulieu et al. 2006), are marked for reference. This 
plot shows that an Earth-mass planet at 1 AU might be detectable if the ring contains grains 
more than a few tens of microns in size and a planet with mass equal to a few times that 
of Mars might be detectable near 10 AU if the ring contains grains more than one hundred 
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microns in size. 

The detectability of a ring structure depends upon the size distribution of dust within 
the ring structure. Dust produced according to a crushing law less steep than the Dohnanyi 
(1969) crushing law (a < 3.4) will result in more highly contrasted ring structures because 
of the increased relative contribution of the large grains. For crushing laws with a = 3 and 
a = 2, the curves of constant maximum particle size shown in Figure 12 shift downward by 
a factor of approximately 1.25 and 1.55, respectively. 

These values are subject to the assumptions of our simulations, which do not include a 
dynamically hot component in the dust cloud. This component would reduce the contrast 
in the ring, making planets harder to detect for a given cloud mass. So the detection limits 
shown in Figure 12 should be thought of as best-case scenarios. 


5. Caveats 

Our simulations include a number of simplifying assumptions, which we summarize here. 
We ignored the effects of dynamically hot dust, like dust that might come from comets. 
Trapping probability decreases dramatically for particles on highly eccentric and inclined 
orbits, so we expect dynamically cold dust to dominate any resonant debris disk structure. 
As a first approximation, we can treat the contribution from the dynamically hot dust as 
a constant surface density cloud component, which reduces the contrast of any structure 
formed from the dynamically cold component. Estimates of the ratio of asteroidal dust to 
cometary dust in our solar system range from 1:10 to 7:10 (Ipatov et al. 2007). For other 
systems, this ratio is also unknown. 

Our simulations also assumed a single planet on a circular orbit around a Sun-like star. 
We have performed trial simulations of our solar system and demonstrated that the presence 
of Jupiter may reduce the Earth’s ring contrast. Other multiple-planet systems may also 
exhibit a similar effect. Additionally, planets on eccentric orbits give rise to additional MMRs 
with different capture probabilities and geometries (Kuchner & Holman 2003). 

The ring contrasts of inclined systems can vary significantly depending on the inclination 
and radial extent of the dust cloud. In edge-on systems, resonant features can overlap as 
seen from the Earth, complicating their interpretation. The contrasts we provide are useful 
only to systems for which projection effects can be taken into account. 

Finally, our multi-particle-size models demonstrate the subtlety of collisional effects in 
dust clouds. Collisional effects can determine the relative populations of large and small 



grains and potentially alter the morphology of the ring structures. Our simulations can not 
yet handle these effects in detail. 


6. Conclusions 

We have implemented our own hybrid symplcctic integrator for the n-body problem and 
used it to simulate collisionless debris disks, taking into account solar wind and drag effects. 

Each simulation contained 5,000 particles. We found that this number of particles suffices 
to populate the dominant MMRs of a low-mass planet with an accuracy at the few percent 
level, yielding for the first time models of the surface brightness distributions of exozodiacal 
clouds that we can use to quantitatively study the contrasts of resonant features — not just 
their geometries. 

We generated a catalog of resonant structures induced by a single planet on a circular or- 
bit around a Sun-like star, available online at http:/ / asd.gsfc.nasa.gov/Christopher.Stark/catalog.php. 
We investigated 120 sets of model parameters, spanning a range of planet masses, planet 
semi-major axes, and values for f3, assuming dust grains launched from orbits with low ed us t 
and idust- The resulting ring structures exhibited leading-trailing asymmetries, gaps near the 
locations of the planets, and sharp inner edges at ~ 0.83a p . 

We performed a detailed analysis of the surface density contrasts of the rings (Figure 
9). We showed that for a planet on a circular orbit, the contrast and morphology of the 
rings are to good approximation functions of only two parameters, M p and yTqj//?, for a 
given stellar mass and distribution of dust sources for M p < 5 M ffi and f3 < 0.25. Equation 4 
summarizes the contrasts of our single-particle-size models as a function of these parameters. 
Considering only the dynamically cold particles analogous to particles released by asteroids 
in the solar system, we find that terrestrial-mass planets are capable of producing resonant 
ring structures with azimuthally averaged contrasts up to ~ 7 : 1. 

By combining our simulations of grains with particular f3 values, we assembled multi- 
particle-size models of 25,000 particles each. Releasing the particles according to a Dohnanyi 
(1969) crushing law without any subsequent collisional processing results in composite clouds 
whose optical depths are dominated by large particles; large particles will dominate images 
of these clouds in visible light and throughout the IR. Based on these composite models, 
we suggested that the best current models for exozodiacal clouds are those with a narrow 
range of grain sizes corresponding to grains whose collision time roughly equals their PR 
time. Future models should account for processes like grain-grain collisions that destroy 
large grains. 
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Equations 14 and 15 provide semi-analytic predictions for the contrast in optical depth 
of a multi-particle-size cloud of dynamically cold grains. For ring structures composed of 
silicate grains released according to a Dohnanyi crushing law (a = 3.4), Equation 15 gives 
an approximate contrast of 
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where (C t; aa,ie) is the ratio of the azimuthally averaged optical depth in the ring structure 
to the azimuthally averaged background optical depth, .s max is the maximum grain size in 
the ring structure, and s 2 is given by Equation 11. For the case s max > s 2 , the first two 
terms in Equation 16 represent the contrast in the adiabatic limit. The remaining term (and 
the terms in the s max < s 2 case) represents deviations from this limit for smaller particles or 
smaller semi-major axes. 


We plotted the mass of the smallest planet that could be detected through observation 
of a resonant ring structure as a function of planet semi-major axis and particle size in 
Figure 12. We assumed a cloud composed of a range of particle sizes adhering to a Dohnanyi 
(1969) crushing law and a minimum detectable optical depth contrast of 1.5:1. We found 
that planets with masses just a fraction of the Earth’s may form detectable ring structures 
if the rings harbor grains more than several tens of microns in size. 
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Fig. 1. — Maximum fractional error in energy during a 3,000-year integration as a function of 
perihelion distance for two scenarios: a two-body system of the Sun & Jupiter (solid line) and 
a three-body system of the Sun, Jupiter & Saturn (dashed line). For the two-body system, 
Jupiter’s perihelion distance is plotted. For the three-body system, Saturn’s perihelion 
distance was altered while Jupiter’s remained fixed. The inclinations and eccentricities of 
both planets remained fixed, cf. DLL98 Fig. 3. 
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Fig. 2. — Top: eccentricity as a function of time for a dust particle with /3 = 0.2 and 
sw = 0.35. Bottom: semi-major axis as a function of time. Our results match those of 
MMM02 (cf. MMM02, Fig. 1) and agree with the analytical solution (Wyatt & Whipple 
1950). 
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Fig. 3. — The integrated trajectory of comet P/Oterma during a close encounter with Jupiter 
as viewed in the frame centered on and rotating with Jupiter with the Sun on the negative 
rc-axis. Shown are the results of a Bulirsch-Stoer integrator and our hybrid symplectic 
integrator for four values of integration time step. The hybrid symplectic results overlap the 
Bulirsch-Stoer results for a timestep of 1 day (cf. C99, Fig. 4). 




25 



(Linear Scale) (Linear Scale) -7.2 -3.2 0.8 4.8 8.8 

Surface Density Surface Density a 


Fig. 4. — Comparison of our hybrid symplectic integrator with a Bulirsch-Stoer integrator. 
Left: surface density histogram for 1,000 particles in the Sun- Earth system using a Bulirsch- 
Stoer integrator. Middle: surface density histrogram for the same initial conditions using our 
hybrid symplectic integrator. Right: Bulirsch-Stoer histogram minus the hybrid symplectic 
histogram (image is in units of a, the \pn Poisson noise associated with the histograms). 
Except in a handful of pixels, the difference is roughly consistent with Poisson noise. 
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Fig. 5. — Population of Neptune’s MMRs for three independent simulations of 5,000 particles 
each (shown in green, red, and black). Populations of the 2:1 and 3:2 resonances differ among 
the three simulations by 6.4% and 4.3%, respectively (c.f. Moro-Martfn & Malhotra 2002, 
Fig. 5). 
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Fig. 6. — Surface density distributions for four of the 120 simulations (scale is relative). 
The star is located at the center of the image and the planet is marked with a white dot. 
The planet orbits counter-clockwise in these images. Integrations were truncated at half the 
planet’s semi-major axis. The simulations shown on the right have different values of a p and 
(3, but the same value of v /ajj//3. Their surface density distributions are nearly identical; 
their difference is consistent with Poisson noise (see Section 3.2). 
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Fig. 7. — Left: the angular size of the “gap” in the ring structure around the location of the 
planet in our simulations versus the contrast of the ring structure. Right: the sine of the 
prograde shift plotted against the function (3(1 — (3)/(M p a p ). We removed all data with 
Caa,ie <1-6 from these plots. Solid lines show linear fits to the data. 
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Fig. 8. — The azimuthally-averaged contrast measured at the inner edge of the ring structure 
(see Section 3.2 for definition of contrast) as a function of (3 (left panel) and planet mass 
(right panel). Both figures show a transition from a no-trapping regime to a saturation 
regime where contrast is independent of semi-major axis. 
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Fig. 9. — The contrast in surface density of the ring structure compared to the background 
cloud for all combinations of M p , a p , and f3 (see Section 3.2 for definitions of contrast). The 
contrast is only a function of two parameters: planet mass and The solid lines are 

fits to the data (see Equation 4). 
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Fig. 10. — Comparison of the optical depths for a composite cloud formed by two different 
methods for M p = 2.0 M® and a p = 6.0 AU. The planet, marked with a white dot, orbits 
counter-clockwise in these images. Top-left: A composite collisionless cloud where the parti- 
cles are released with a size distribution equal to the crushing law used by Dohnanyi (1969). 
Top-right: The same composite cloud, but formed by forcing the surface density outside of 
the ring structure to obey a Dohnanyi distribution. Bottom-left: The optical depth of the 
smallest particles included in the composite clouds. Bottom-right: The optical depth of the 
largest particles included in the composite clouds. The largest particles dominate the optical 
depth in a cloud of particles released with a Dohnanyi crushing law. 
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Fig. 11. — Contrast in optical depth for multi-particle-size clouds (solid lines) compared to 
single-particle-size clouds (dashed lines) assuming a Dohnanyi (1969) crushing law (a = 3.4; 
see Equations 16). From top to bottom, the six solid lines and six dashed lines correspond 
to six values of planet mass: 5, 2, 1, 0.25, and 0.1 M e . The contributions of the small grains 
reduce the contrasts of the multi-particle-size clouds compared to single-particle-size clouds 
with the same minimum value of (3. 





Fig. 12. — Minimum detectable planet mass in a multi-particle-size collisionless cloud as a 
function of semi-major axis and maximum grain size, assuming a Sun-like star, a minimum 
detectable ring contrast of C ri AA,iE = 1.5, and dust produced according to a Dohnanyi (1969) 
crushing law (a = 3.4). Earth-like and Mars-like planets are denoted with an E and M, 
respectively. The 5.5M ffi exoplanet OGLE-2005-BLG-390Lb is denoted with an O (Beaulieu 
et al. 2006). Listed values for maximum dust size in the ring structure assume perfectly 
absorbing spherical grains with mean density p = 2.0 gm cm -3 and radius s max . The bold 
line shows the case of the solar zodiacal cloud, for which the observed emission is dominated 
by 30 pm grains (Fixsen & Dwek 2002). The dashed lines show typical inner and outer 
detection limits for a mission similar to TPF. 




